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ABSTRACT 


The  ability  to  measure  and  compensate  for  power  line  harmonics  has  become  a 
growing  area  of  interest  because  of  today’s  commonly  used  electronic  equipment.  Since 
the  number  and  relative  magnitudes  of  the  harmonics  on  the  power  line  are  a  function 
of  the  load,  estimation  of  an  equivalent  load  can  be  accomplished.  Because  of  variation 
in  the  load,  the  need  for  an  adaptive  algorithm  is  imperative.  Thus  far,  few  algorithms 
for  determining  harmonic  contents  have  not  dealt  with  the  problem  associated  with  the 
limited  power  of  the  line  conditioner. 

This  thesis  deals  with  a  previously  known  harmonic  compensating  algorithm  and 
introduces  a  new  algorithm  which  both  compensates  for  harmonic  noise  and  estimates  the 
load  as  a  transformation  matrix  depending  on  the  associated  transfer  function  of  the 
active  line  conditioner  in  use. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Many  of  today’s  electronic  devices  /. e. ,  computers,  fluorescent  lights  and 
microwave  ovens,  effect  power  distribution  due  to  their  nonlinear  consumption  of  power. 
The  result  is  irregular  current  and  voltage  wave  forms  on  the  power  line.  Active  line 
conditioners  provide  a  way  of  eliminating  the  accompanying  noise  on  a  power  line  by 
independently  adjusting  the  active  and  nonactive  components,  thereby  maintaining  a 
constant  sinusoidal  bus  voltage. 


Figure  1.1  Bus  with  Linear  and  Nonlinear  Loads  and  Active  Line  Conditioner 

The  amount  of  voltage  distortion  is  a  function  of  the  nonlinear  load  distribution, 
their  current  spectra,  topology  and  frequency  dependance  within  the  network.  For  a 
nonsinusoidal  voltage  [Ref  1] 

VB  =  y^KjSincof  +  v/2£  ^sin (Auf  +  <xh)  (1.1) 

h*l 


1 


with  linear  and  non  linear  loads  drawing  a  total  current. 


where 


h  = 


*al 


+  lrI  +  l 


a 


ial  =  y/5/iCOS01sinQr 

is  the  active  component  of  the  fundamental  current, 


(1.2) 


(1.3) 


irl  =  v/2i1sin01coso>r 


(1.4) 


is  the  reactive  component  of  the  fundamental  current  and 


ig  =  v/25^/*sin(Atof  +  ah  +  dA) 

h*  1 


is  the  harmonic  current.  The  apparent  power  can  be  described  as  follows. 


(1.5) 


S  =  V  I 

rmsrms 


N 


W  *  ♦  E7*) 

A#1  A*1 


(1-6) 


Equation  1.6  can  be  expressed  in  phasor  form, 


s  -  M  *  P*)2  ♦  Q?  ♦  Q 


.2 

Jf 


(1.7) 


where 

P,  =  VlIlcosQh  (1.8) 

represents  the  Fundamental  Power  Frequency  Active  Power.  The  Harmonic  Active 
Power  is, 
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(1.9) 


<■„  -  E 

AM 


and 


Q,  =  Kj/jSmOj  (1.10) 

is  the  Power  Frequency  Reactive  Power.  The  Harmonic  Reactive  Power  is 

QH  =  HVHIk^nQH  (1.11) 

A*  1 


The  components  of  Qw  are  generated  by  specific  harmonic  voltages  and  harmonic 
currents  (in  no  particular  order)  [Ref  2].  Because  of  the  vector  properties  of  the  reactive 
power  components,  control  or  cancellation  is  possible  using  vectors  of  identical  frequency 
and  magnitude  with  opposite  phase.  Therefore,  the  Harmonic  Reactive  Power  can  be 
eliminated  by  introducing  or  drawing  current  from  the  power  line  which  is  180°  out  of 
phase  with  each  respective  harmonic.  With  this  in  mind,  the  conditioner  current  is  as 
follows, 


where  i^,  and  iCrl 
respectively,  and 


lC  =  lCal  +  lCrl  +  lCB 


(1.12) 


are  the  conditioner  Fundamental  Active  and  Reactive  currents 


iCH  =  4E /cAsin(/rwf  +  «a  +  Y  H)  (1.13) 

A-l 

is  the  conditioner  harmonic  current. 


3 


B.  ACTIVE  LINE  CONDITIONERS 

Active  line  conditioners  serve  a  dual  purpose.  First,  they  adjust  one  or  more  loads 
thereby  changing  the  active  power.  Secondly,  they  are  capable  of  controlling  the 
amplitude  and  phase  characteristics  of  the  nonreactive  currents  iCrl  and  iCH  which  affect 
the  value  of  Q,  and/or  QH  [Ref  2,4].  By  using  a  solid  state  switching  network  at  a 
frequency  much  greater  then  that  of  the  fundamental,  the  line  current  is  modulated  in 
order  to  maintain  the  boundaries  of  a  desired  template  z,  shown  in  Figure  1.2.  For  a 
narrow  boundary  of  error  5,  the  conditioner  current  i'c  is  unaffected  by  the  fluctuations 
(the  zig-zaging)  within  the  boundary  area.  By  adjusting  the  template  waveform  through 
a  feedback  circuit,  the  conditioner  current  spectrum  can  be  altered  to  produce  an  overall 
current  tr  that  is  as  close  to  sinusoidal  as  possible.  To  do  this,  the  load  current  iL  is 
monitored  to  determine  the  harmonics  present.  Then  by  injecting  a  current  ic  from  the 
line  conditioner  which  is  180°  out  of  phase  with  that  of  the  load,  the  unwanted  harmonics 
can  be  cancelled  out. 


1.  Equivalent  Circuit  Modeling 

For  medium  and  low  voltage  systems,  the  best  practical  means  of  adjusting 
the  conditioner  current  ic  is  by  minimization  of  the  voltage  distortions  at  the  conditioner 
location  on  the  bus  [Ref.  5].  The  voltage  at  the  conditioner  node  can  be  represented  by 


<‘-w> 

1-1 

where 

1^  =  Bus  x  harmonic  current  phasor  of  order  h 


4 


template  z 


Figure  1.2  Bus  Voltage,  Line  Conditioner  Template  and  Current 

=  Harmonic  Complex  Impedances  of  the  node 
N  =  Number  of  independent  nodes 

A  Norton  equivalent  circuit  for  the  bus  has  the  equivalent  harmonic  current 

I.»  -  E  ZBxklxX  (115) 

x-1  x*B 

where  YBBh  =  l/ZBBh  is  the  self-admittance  of  the  node  for  the  harmonic  h.  Figure  1.3 
shows  the  Norton  equivalent  circuit  with  the  associated  harmonic  currents  and  load  Za. 


Figure  1.3  System  Approximation  using  Norton  Equivalent  Circuit 
Where  is  given  in  Equation  1.16, 
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^sa  ~  &sh  *  jk<*LSh  -  ZBaAJZaA 


(1.16) 


and  %bh  is  the  load  impedance  of  the  harmonic  order  h. 

From  the  equivalent  load  Zsa  the  voltage  due  to  the  offending  harmonics  \Hh 
can  be  defined  as 

V**  - 1 A  (1.17) 

For  a  linear  resistance  and  inductor  and  L^,,  a  conditioner  current  ic  equal 
to  the  negative  of  IfA  would  eliminate  the  harmonic  voltage  vWA.  It  is  important  to  note 
that  active  line  conditioners  are  inherently  limited  in  their  maximum  current  output, 
therefore,  negating  the  entire  value  of  Irt  may  not  be  possible.  Although  limited,  any 
reduction  of  the  harmonic  noise,  especially  of  lower  order,  significantly  improves  the 
recognition  of  the  fundamental. 
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II.  SURVEY  OF  PERVIOUS  WORK 


A.  ADAPTIVE  ESTIMATION  OF  HARMONIC  VOLTAGE 

The  best  fitting  sinusoidal  wave  to  a  nonsinusoidal  periodic  wave  is  the  fundamental 
[Ref.  1].  The  error  associated  with  such  a  system  can  be  written. 

*  =  vb  ~  vi  =  vh  +  vc  (21> 

Since  the  signals  vH  and  v,  are  periodic,  the  error,  e,  is  statistically  stationary. 
Therefore,  the  expected  value  of  the  square  of  the  error  e  results  in  a  quadratic  function 
which  has  a  guaranteed  minimum  for  real  physical  signals  [Ref.  6].  Then  by  minimizing 
the  mean  square  of  the  error  (MSE),  the  signal  should  be  nearly  identical  to  that  of  the 
fundamental.  By  representing  the  error  voltage  due  to  the  harmonics  as  the  sum  of 
weighted  sines  and  cosines,  an  error  surface  for  each  weight  can  be  defined  thereby 
making  the  calculation  of  a  minimum  possible.  Figure  2.1  shows  the  block  diagram  for 
such  an  adaptive  system. 

*S2,  *S3  •"  *sh  and  xa>  x&  -  *ch  represent  discrete  versions  of  the  harmonics 
associated  with  the  power  line,  and  C2,  C3,  -  Ch  and  S2,  S3,  Sh  are  the  weights  of  the 
associated  harmonic.  The  correction  voltage  to  the  conditioner  is  defined  by 

yk  =  [Shsm(ho)kT/N)  +  Chcos(h(okT/N)]  (2.2) 

h*  1 

where  k  =  time  index 

T  =  l/f  =  2-k!<j3  =  period  or  one  cycle 
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Figure  2.1  Block  Diagram  of  Adaptive  System 

N  =  number  of  samples  per  cycle 
The  Active  Line  Conditioner  current  is  given  by 

iCH  =  KD^  [5Asin(A(of)  +  Chcos(h(ot)] 

h*  1 

where  D  =  The  gain  of  the  D/A  converter 

K  =  Converter  constant  of  the  Conditioner  in  (A/V) 

The  discrete  error  ek  as  a  function  of  voltage  becomes 

**  =  VC*  +  VHk 

where  va  and  vm  are  represented  by 

vc*  =  KDH  ^Sh[Shs\n{hu)kTIN)  +  Chcos(hi»kTIN)] 

h*  1 
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(2.6) 


vHk  =  v^2^  ZAI<A[cospAsin(/i(oik7/^)+  sw$hcos(h<jjkT/N)] 

h*l 

Simply  setting  the  sine  and  cosine  weights  equal  to 


s„  -  -v2(/;,cosp,)//fD 
c,  =  -v5(/rtcosp,)/A:D 


requires  the  error  to  be  zero.  Since  the  actual  load  impedance  ZSh  is  not  known  and 
changes  with  time,  an  estimation  of  the  sine  and  cosine  weights  is  performed  by  the  MSE 
processor  using  the  following  linear  prediction. 

^*+1  =  +  (~Ys )n/^  (2  8) 

Ck.i  =  M-Vc)p//i 

where 

Ys  =  5e/3S  /2  q\ 

vc  =  de/ac  c  ■ 

are  the  error  gradients  of  the  sine  and  cosine  weights  respectively,  and  n  is  a  constant 
called  the  acceleration  factor  which  is  directly  related  to  both  the  rate  of  convergence, 
and  the  magnitude  of  any  over-shoot  in  reaching  the  minimum.  The  h  term  is  used  to 
scale  the  acceleration  factor  by  an  amount  proportional  to  the  harmonic  being  evaluated. 
This  allows  for  faster  convergence  of  lower  order  harmonics,  the  largest  error,  without 
driving  the  higher  order  unstable.  Figure  2.2  shows  a  quadratic  error  surface  as  a 
function  of  a  single  weight. 


9 


- 1 - 1 - 1 - I - 

Ck  Ck+1  ck+4  Ck+3 

WEIGHT 

Figure  2.2  Single  Weight  Error  Surface  Example 

Starting  on  the  left  side  where  Ck+I  >  Ck  and  ek+1  <  e*  the  gradient  is  negative 
and  the  value  of  e  converges  toward  the  minimum.  If  the  value  of  the  weighting  factor 
produces  a  higher  error  then  the  previous  £k+<  >  ek+i,  an  over  shoot  occurs,  but  the 
gradient  remains  negative  thereby  predicting  a  smaller  weighting  value  then  the  current 
one  and  £  again  converges  toward  the  minimum.  The  same  result  would  be  obtained  for 
an  initial  weight  greater  then  that  need  to  minimize  e. 

B.  LINEAR  LOAD  SIMULATION 

A  program  for  testing  the  validity  of  the  adaptive  algorithm  was  written  for 
MATLAB  using  the  parameters  in  Figure  2.3  and  a  noise  component  equating  to  ten 
percent  Total  Harmonic  Distortion  (THD)  assuming  only  odd  harmonics  up  to  the 
twenty-first  [Ref.  1].  The  fundamental  frequency  is  assumed  to  be  a  standard  60  Hz 
outlet  and  the  line  conditioner  to  have  no  power  limitation. 
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Figure  2.4  shows  the  MSE  for  the  first  three  offending  harmonics  and  their 
convergence  to  zero,  all  well  under  one  second  using  an  acceleration  constant  of  3e-7. 
The  THD  in  Figure  2.5  also  follows  a  similar  pattern  since  it  is  most  affected  by  the 
lower  order  harmonics  and  becomes  one-one  hundredth  of  its  original  value  after  just  one 
half  second,  or  thirty  iterations.  Figure  2.6  shows  the  weighting  coefficients  of  the  sine 
and  cosine  for  the  first  three  offending  harmonics.  With  the  help  of  a  simple 
trigonometric  identity,  these  values  can  easily  be  shown  to  correspond  to  the  magnitude 
and  phase  of  the  original  noise  components.  Also  note  that  the  final  weighting  values 
where  reached  asymptotically  without  any  over-shoot  indicating  the  choice  of  the 
acceleration  coefficient  is  optimal  with  respect  to  requiring  minimum  power  from  the 
active  line  conditioner. 
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Magnitude  (volts) 


Mean  Square  Error 
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Total  Harmonic  Distortion 
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Magnitude 


Time  (sec) 


Figure  2.6  Sine  &  Cosine  Weighting  Coefficients  for  Linear  Load 
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C.  NONLINEAR  LOAD  SIMULATION 


The  inductor  in  Figure  2.3  was  substituted  for  the  one  shown  in  Figure  2.7  to 
produce  a  nonlinear  response  in  the  load.  The  nonlinearity  was  chosen  to  provide 
approximately  ten  percent  deviation  from  that  of  the  linear  case. 


The  nonlinear  load  results  were  very  similar  to  those  from  the  linear  with  small 
deviations  in  the  THD  and  identical  results  for  the  weighting  coefficients.  Some  of  these 
values  are  summarized  below  in  Table  2.1. 

TABLE  2.1  TOTAL  HARMONIC  DISTORTION  FOR  LINEAR  AND  NONLINEAR 
LOADS 


Time  (sec) 

.0167 

.0333 

.0667 

.133 

.250 

.500 

1.00  1 

Linear 

9.327 

8.849 

6.381 

3.300 

1.047 

.123 

.0011 

Nonlinear 

9.245 

8.874 

6.410 

3.427 

1.120 

.120 

.0012 
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Mean  Square  Error 
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Total  Harmonic  Distortion 
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Magnitude 


Weighting  Coefficients 


Figure  2.10  Sine  &  Cosine  Weighting  Coefficients  for  Nonlinear  Load 
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ffl.  OPTIMAL  ESTIMATION  WITH  SYSTEM  IDENTIFICATION 


A.  MODELING  THE  NETWORK 

The  model  used  for  optimal  estimation  is  very  similar  to  that  found  in  Figure  2.1 
with  the  exception  that  the  impedance,  although  unknown,  will  be  estimated  along  with 
minimizing  the  harmonic  error.  The  block  diagram  model  is  shown  in  Figure  3.1  with 
the  impedance  of  the  present  load  represented  as  a  square  matrix  H.  It  should  be  noted 
that  this  model  can  be  used  for  different  system  transfer  function  input/output 
relationships  other  than  current  to  voltage. 


w(n) 


Figure  3.1  Optimal  Estimation  Impedance  Model 

Since  the  harmonic  noise  w(t)  is  assumed  to  be  harmonics  of  the  fundamental  it  can 
be  represented  as  a  sum  of  sinusoids  shown  in  Equation  3.1  for  a  continuous  signal. 


(3.1) 


" k 

w(t)  =  £  Axcos(2icx/O0  +  Bxsm(2xx/0f) 

x-2 


For  a  discrete  sampled  system  Nyquist  criteria  must  be  maintained,  therefore  the 
sampling  frequency,^,  must  be  an  integer  value  of  the  fundamental  which  is  greater  then 
twice  the  highest  harmonic  frequency  to  be  eliminated.  The  continuous  frequencies  are 
converted  to  discrete  under  the  following  conditions. 

2*/—  =2 it—  where  1  z  l  z  Nh  <  —  (3.2) 

fs  M  *  2 

From  this  the  disturbance  w(t)  can  be  represented  discretely  in  matrix  form  as 

w(/i)  *  WTx(n)  0-3) 

where 


WT  -  [0,  0,  Bji  Bj,  ■  ■ 

112  2 
x(#i)  =  [cos2it — n,  sm2x  —  n,  cos2ic  —  n,  sin2ii — n,  •  •  •  , 
M  M  M  M 

Nk  Nh  . 

cos2x  —  n,  sm2rt — n] 
M  M 


(3.4) 


Because  the  harmonic  noise  does  not  change  instantaneously  with  changes  to  the 
load,  it  is  reasonable  to  assume  that  it  is  periodic.  By  dividing  the  time  scale  into 
periodic  intervals  of  length  N,  which  is  a  multiple  of  M,  then  for  all  n  x(n)  =  x(n  +  M) 


t 


N  (k  +  1)N  (k  +  2)N 

Figure  3.2  Discrete  Time  Period 

By  defining  the  control  input  current  to  be 

u(n)  =  u*x(n)  jfcNsusitN  +  N-1  (3-5). 

where  u*  is  a  weight  vector  to  be  determined.  The  error  can  now  be  written  discretely 
in  terms  of  x(n)  as 

e(n)  =  WTx(n)  +  ujHTx(n)  (36) 

Note  that  for  a  linear  impedance  H  becomes  a  diagonal  matrix,  otherwise  it  is  not. 

Since  w (n)  is  the  signal  which  is  to  be  eliminated,  by  using  its  frequency 
information  in  the  control  input,  u(n),  it  will  provide  some  of  the  needed  information  to 
negate  it.  The  remaining  control  information  will  come  from  a  recursive  estimation  of 
its  Fourier  coefficients  which  is  the  main  basis  of  this  algorithm. 

B.  CONTROL 

Since  the  value  of  H  can  be  estimated  through  the  use  of  system  identification 
techniques,  it  can  be  assumed  to  be  known  for  the  purpose  of  determining  the  control 
necessary  to  eliminate  the  harmonic  noise.  From  Equation  3.6  it  is  easy  to  see  that  WT 
and  utTHT  are  scalars  which  can  be  combined  to  represent  the  error  weights  for  each 
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sample  point  over  a  single  interval,  N,  of  the  fundamental.  The  error  associated  with 
a  single  sample  n  is  shown  in  Equation  3.7. 


e(n)  =  ejx(n)  0-7) 

Now  each  term  in  Equation  3.6  can  be  defined  in  terms  of  x(n). 

Due  to  the  periodicity  of  the  harmonic  frequencies,  the  frequency  components  can 
be  eliminated  by  subtracting  the  error  of  the  k*  interval  from  the  (fc-l),h  interval  resulting 
in  a  difference  equation  of  just  weighted  vectors  shown  below. 

e*x(/i)  =  WTx(n)  +  u*HTx(n ) 

eJ.,x(n-N)  =  WTx(/i -N)  +  utT.,HTx(w-N)  (3.8) 

ej  -  «£,  =  HT(nJ  -  «£,) 

After  some  manipulation,  Equation  3.8  can  be  written  as 

Q(e*  ~  Vi>  =  '  u*-i  <3-9> 

where  Q  =  H'1 

Now  that  the  control  is  in  terms  of  the  error  and  an  admittance  matrix  Q,  using  a 
linear  predictor  similar  to  that  in  Equation  2.8  can  be  used. 

=  u*-i  -  «Qe*-i  <3-10> 

where  a  is  a  scalar  defmed  on  the  interval  -1  <  a  <1. 

The  error,  e*,  can  be  driven  to  zero  for  Q  not  equal  to  the  null  set. 
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C.  ESTIMATION 


Because  the  load  to  the  bus  changes  over  time  some  method  of  updating  the  value 
of  H  to  those  changes  must  exist  in  order  for  the  controller  to  effectively  reduce  the 
harmonic  noise  present.  Estimation  of  the  admittance  matrix  Q  can  easily  be 
incorporated  with  the  control  through  system  identification  methods  using  a  recursive 
least  squares  (RLS)  algorithm.  By  choosing  estimates  of  the  control  and  error  to  be 

“*  =  U*  '  Utl  (3.11) 

-  e*  "  e*-i 

Equation  3.9  reduces  to  a  matrix  form  of  the  RLS  equation. 

u,  -  ejQ  a  12) 

Although  Equation  3.12  is  in  matrix  form,  it  is  important  to  remember  that  the 
estimate  of  each  row  of  Q  is  a  unique  difference  equation  of  the  associated  control  and 
error  coefficients  of  their  respective  frequency.  In  other  words,  even  though  the 
frequency  vector  x(n)  is  not  found  in  Equation  3.12  the  relationship  between  the  third 
harmonic  error  and  control  coefficients  remains  linear.  It  is  well  known  that  the  output 
of  a  linear  system  differs  from  the  input  only  in  magnitude  and  phase,  therefore  the 
system  output  y(t)  can  be  written  in  terms  of  u(t)  as  follows: 


yk(0  =  |Htj  Akcos(2it*/ar  +  at)  +  |Ht|Bksin(27tfc/or +  at)  = 

Ckcos(27tfc/or)  +  Dksin(27t  fc/0f) 


(3.13) 


where  Ck  and  reflect  the  magnitude  and  phase  changes  of  the  system  on  the  input 
coefficients  A*  and  B*.  With  the  help  of  a  trigonometric  identity.  Equation  3.13  can  be 
written  in  a  more  convenient  matrix  form  [Ref  7). 


cosak 

sinat 


-sinat 

cosak 


(3.14) 


Since  H*,  cosa*  and  sina*  are  scalars  Equation  3.14  can  be  arrange  in  a  recursive  form 
similar  to  that  of  Equation  3.12  as  follows 


Ck 

V 

Bk 

Yk 

Dk 

Bk 

Ak 

(3.15) 


where  Yk  and  Zk  represent  the  estimate  for  the  product  of  the  magnitude  and  phase 
change  for  the  cosine  and  sine  terms  respectively  of  a  given  harmonic.  Equation  3.15 
represents  the  estimates  of  the  coefficients  for  just  a  single  harmonic  of  the  system  and 
can  be  thought  of  as  a  building  block  of  the  matrix  RLS  Equation  3.12. 

Now  all  the  needed  information  is  available  for  implementation  of  a  Kalman  filter 
based  RLS  estimation  [Ref  8,  9].  Figure  3.3  shows  a  block  diagram  representation  of 
the  system  model  with  the  estimation  algorithm. 

1.  Considerations  in  Applying  FFT  to  Harmonic  Analysis 

Since  this  algorithm  emphasizes  the  use  of  Fourier  coefficients,  some  aspects 
of  using  the  Fast  Fourier  Transform  (FFT)  will  be  addressed.  The  FFT  algorithm  has 
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w 


Figure  3.3  System  Diagram  with  Estimation  and  Control 

useful  applications  in  power  system  networks,  but  can  produce  erroneous  information  if 
not  applied  correctly.  Certain  assumptions  about  the  FFT  must  be  understood  to  avoid 
false  representation  of  the  associated  signal  [Ref  7]. 

•  The  signal  is  stationary  (constant  magnitude). 

•  Each  frequency  in  the  signal  is  an  integer  multiple  of  the  fundamental. 

•  The  sampling  frequency  is  equal  to  the  number  of  samples  times  the  fundamental 
frequency  used  in  the  algorithm. 

•  The  sampling  frequency  meets  Nyquist  criteria. 

D.  SIMULATION  RESULTS 

In  testing  the  optimal  estimation  algorithm  the  same  harmonic  noise  components 
from  the  MSE  in  Chapter  II  were  used.  Since  the  impedance  matrix  of  the  system  is 
estimated  using  RLS  any  linear  transformation  matrix  for  H  can  be  used.  For  the 
purpose  of  this  research  a  simple  diagonal  matrix  with  a  linear  progression  from  one  to 
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forty-two  was  used.  Figure  3.4  shows  the  control  input  to  the  line  conditioner  for  the 
three  highest  offending  harmonics.  As  in  the  MSE  case,  the  values  are  reached 
asymptotically  without  any  overshoot,  thus  showing  the  stability  of  the  algorithm.  The 
optimal  estimation  algorithm  demonstrated  superior  robustness  and  stability  compared  to 
that  of  the  MSE  with  respect  to  the  linear  predictor  constant  o.  The  optimal  estimator 
provided  stable  and  consistent  results  for  positive  values  of  a  up  to  one.  While 
individual  harmonics  in  the  MSE  case  were  highly  sensitive,  and  often  grew  unstable, 
with  changes  in  a. 
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Control  Weighting  Coefficients 


Time  (sec) 


Figure  3.4  Control  Input  to  Conditioner  U* 
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Figure  3.5  shows  the  total  harmonic  distortion  for  several  values  of  a  as  a  function 


of  time. 


>i  v  ^ _ i _ i  .  „  —  _  .1,1 

0  0.2  0.4  0.6  0.8 


Time  (sec) 

Figure  3.5  Total  Harmonic  Distortion  (THD) 
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An  additional  piece  of  information  provided  by  the  optimal  estimation  algorithm 
is  an  adaptive  estimate  of  the  load  impedance  in  the  form  of  a  matrix,  H.  Table  3.1 
gives  a  break  down  of  several  of  the  actual  and  estimated  matrix  values  along  with  their 


respective  errors. 

TABLE  3.1  ACTUAL  AND  ESTIMATED  IMPEDANCE  MATRIX  COEFFICIENTS 


Harmonic 

3rd 

5th 

7th 

9th 

11th 

Actual 

6 

a 

10 

11 

14 

15 

18 

19 

22 

23 

Estimated 

6 

a 

10 

11 

14 

15 

18 

19 

22 

23 

Percent 

error 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

It  is  important  to  point  out  that  only  harmonic  frequencies  which  are  present  in  the 
offending  noise  will  produce  impedance  coefficient  estimates  for  the  H  matrix.  This  is 
simply  the  result  of  having  no  error  to  drive  the  RLS  equation.  If  noise  other  than  an 
odd  harmonic  of  the  fundamental  were  present  an  estimate  of  that  impedance  would 
appear  in  the  matrix  H. 


rv.  CONCLUSIONS 


A.  MSE  ALGORITHM  VS.  OPTIMAL  ESTIMATOR 

When  the  load  is  linear,  the  optimal  estimation  algorithm  proved  to  be  much  more 
effective  in  eliminating  the  associated  harmonic  noise  and  more  robust  with  respect  to 
changes  in  the  gain  of  the  linear  prediction,  as  indicated  in  Equations  2.8  and  3. 10.  In 
addition  to  elimination  of  the  harmonic  noise,  the  optimal  estimator  provides  a  linear 
representation  of  the  present  system  load  in  the  form  the  impedance  matrix  H.  The 
impedance  information  of  the  H  matrix  is  beneficial  in  helping  to  determine  the  proper 
specifications  of  the  active  power  line  conditioner  for  the  particular  application. 

B.  FUTURE  RESEARCH 

Because  the  optimal  estimation  algorithm  uses  the  Fast  Fourier  Transform  of  the 
error  signal  it  is  currently  limited  to  the  linear  case.  Nonlinearities  in  the  error  signal 
present  a  difficult  obstacle  to  overcome  using  standard  Fouier  transforms.  Investigation 
into  adapting  the  optimal  estimation  algorithm  for  nonlinear  contingencies  would  be 
highly  beneficial  and  would  provide  a  better  and  more  accurate  model  of  the  power  line 
load  impedance. 
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APPENDIX  A  -  (MATLAB)  ADAPTIVE  LINEAR  MODEL 


%  %  THESIS  PROGRAM  I 

%%  JOEL  ZUPFER 

%%  17  MAY  93 

%  %  REVISED  18  NOVEMBER  93 

%  %  SIMULATION  OF  CIRCUIT  FIGURE  2.3 


clear 

R 

L 

f 

N 

lag 

Cycle 

check 

totN 

Delay 

accel 

Wt 

H 

P 

mu 


=  700; 

=  1.857; 

=  60; 

=  120; 

-  120; 

=  input(’Number  of  cycles  =  ’); 

=  inputfNumber  of  cycles  between 

=  N*round(Cycle/2)*check; 

=  N*check; 

=  3e-7; 

=  [.51  .16  .056  .035  .025  .025  .02 


%MODELLED  RESISTOR  VALUE  (OHMS) 
%MODELLED  INDUCTOR  VALUE  (HENERY’S) 
%  FUNDAMENTAL  FREQUENCY  (HZ) 
%  NUMBER  OF  SAMPLES  PER  CYCLE 
%LAG  FOR  CALCULATION  OF  THE  MSE 
96NUMBER  OF  CYCLES  IN  SIMULATION 
gradient  calculations  =  ’); 

%  CYCLES  BETWEEN  GRADIENT  CALCULATIONS 
96TOTAL  NUMBER  OF  POINTS  IN  SIMULATION 
%  COUNTER  FOR  CHECKING  ERROR  VOLTAGE 
%  ACCELERATION  FACTOR 

.02  .015  .015]’; 


%  INITIAL  WEIGHTS  OF  ODD  HARMONICS 
*  [3  5  7  9  11  13  15  17  19  21]’;  %ODD  HARMONICS  OF  FUNDAMENTAL  FREQUENCY 
=  pi*[l/3  1/4  1/5  3/2  4/2  5/2  1/6  2/6  3/6  8/6]’;  %INITIAL  PHASE  OF  ODD  HARMONICS 
=  H.A(-l).*accel;  %  ACCELERATION  FACTOR  ADJUSTED  FOR  HARMONICS 


%%%%%%%%%%%%% %  INITIALIZE  PARAMETER  MEMORY%%%%%%%%%%%%%%%%% 


Ih 

=  zeros(length(Wt),totN); 

Ic 

=  zeros(length(Wt),totN); 

If 

=  zeros(l,totN); 

Ve 

=  zeros(length(Wt),totN-l); 

Vf 

=  zeros(l.totN-l); 

Sh 

=  zeros(length(Wt),round(CycIe/2) +2); 

Ch 

=  zeros(length(Wt),round(Cycle/2)+2); 

MSE 

=  zeros(length(W  t)  .Cycle  + 1); 

thd 

=  zeros(l,Cycle+ 1); 

GradSh 

=  zeros(length(Wt),round(Cycle/2)); 

GradCh 

=  zeros(length(Wt),round(Cycle,'2)); 

%  HARMONIC  CURRENT 
%  CONDITIONER  CURRENT 
%  FUNDAMENTAL  CURRENT 
%  ERROR  VOLTAGE 
%  FUNDAMENTAL  VOLTAGE 
%  HARMONIC  SINE  WEIGHTS 
%HARMONIC  COSINE  WEIGHTS 
96MEAN  SQUARE  ERROR 
%TOTAL  HARMONIC  DISTORTION 
%  GRADIENT  OF  SINE  HARMONICS 
%GRADIENT  OF  COSINE  HARMONICS 


%%%%%%%%%%%%%%%%%%  %INmALCONDmONS%  %%%%%%%%%%»%%%%%%%% 
Sh(:,2)  =  ones(length(Wt), l)./H.*.l;  %INITIAL  WEIGHTING  FACTOR  FOR  SINE  HARMONIC 
Cb(:,2)  =  ones(length(Wt),  1 ).  /H.  *.  1  %  INITIAL  WEIGHTING  FACTOR  FOR  COSINE  HARMONIC 
If(l)  =  10*sqrt(2); 

If(2)  =  9.9*sqrt(2);  %PEAK  CURRENT  VALUES 
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%%%%%%%%%%%%%%%%%% %MAIN  PROGRAM %%%%%%%%%%%%%%%%%%%%%% 
%  %  %  %  %  %  %  %  %  %  %  %BUS  VOLTAGE  WITH  NO  CONDITIONER  CURRENT%  %%%%%%%%%% 
fork  =  3:N-1 

sample  =  H.*2*pi*k/N;  %  DISCRETE  SAMPLE  POINT 

If(k)  =  10*sqrt(2)*cos(2*pi*k/N);  %FUNDAMENTAL  CURRENT 

Ih(:,k)  =  Wt.*cos(sample  +  P);  %HARMONIC  CURRENT  OF  THE  LOAD 

Ic(:,k)  =  <Sh(:,l).*sin(sample)  +  Ch(:,l).*cos(sample));  %CONDITIONER  CURRENT 

Vf(k-1)  =  (If(k)  -  lf(k-2))*(L*N*f/2)  +  If(k-1)*R;  ^FUNDAMENTAL  VOLTAGE 

Ve(:,k-1)  =  (Ih(:,k)  +  Ic(:,k)  -  Ih(:,k-2)  -  Ic(:,k-2)).*(L*N*f/2)  +  (Ih(:,k-1)  +  Ic(:,k-1)).*R; 

%BUS  ERROR  VOLTAGE 
end 

MSE(:,1)  =  sqrt(mean( Ve( : ,  1 :  k- 1 )  ’ .  A2))  ’ ;  %MEAN  SQUARE  OF  BUS  ERROR  VOLTAGE 

thd(l)  =  sqrt(sum(max(V“(:,2:k-l)’).^2/2)’>100/(max(Vf(2:k-l))/sqrt(2)); 

%  TOTAL  HARMONIC  DISTORTION 
final  =  k  + 1 ;  %STEP  INDEX 

%%%%%%%%%%% %%BUS  VOLTAGE  WITH  CONDITIONER  CURRENT% %%%%%%%%%%% 
for  index  =  1  :round(Cycle/2) 

for  k  =  final: Delay  +  final- 1 

sample  =  H.*2*pi*k/N;  %  DISCRETE  SAMPLE  POINT 

If(k)  =  10*sqrt(2)*cos(2*pi*k/N);  %  FUNDAMENTAL  CURRENT 

Ih(:,k)  =  Wt.*cos(sample  +  P);  %HARMONIC  CURRENT  OF  THE  LOAD 

Ic(:,k)  =  (Sh(:,index+  l).*sin(sample)  +  Ch(:,  index).  *cos(sample)); 

^CONDITIONER  CURRENT 

Vf(k-1)  =  (If(k)  -  If(k-2))*(L*N*f/2)  +  lf(k-l)*R;  %FUNDAMENTAL  VOLTAGE 

Ve(:,k-1)  -  (Ih(:,k)  +  Ic(:,k)  -  Ih(:,k-2)  -  Ic(:,k-2)).*(L*N*f/2)  +  (Ih(:,k-1)  +  Ic(:,k-1)).*R; 

%BUS  ERROR  VOLTAGE 
end 


MSE(:,2*index) 

thd(2*index) 

final 

GradShf:,  index) 
Sh(:,index+2) 


=  sqrt(mean(Ve(:,k-l-lag:k-l)’.*2))’;  %MEAN  SQUARE  OF  ERROR  VOLTAGE 
=  sqrt(sum(max(  Ve( : ,  k-  1  -lag  :k- 1)’).^2/2)*1 00/ (max(  Vf(k- 1  -lag  :k- 1  ))/sqrt(2)) ; 

%TOTAL  HARMONIC  DISTORTION 
=  k+1;  %  STEP  INDEX 

=  (MSE(:,2*index)  -  M SE( : , 2 *index- 1 )) . /(Sh( : , index  + 1 )  -  Sh(:, index)); 

%SIN  HARMONICS  GRADIENT  CALCULATION 
=  Sh( index  +  1)  -  GradShf:, index). *mu.*MSE(:,2*index); 

%  PREDICTED  SINE  WEIGHTING  FACTOR 


for  k  =  final: Delay  +  final- 1 

sample  =  H.*2*pi*k/N;  %DISCRETE  SAMPLE  POINT 

If(k)  =  10*sqrt(2)*cos(2*pi*k/N);  %  FUNDAMENTAL  CURRENT 

Ih(:,k)  =  Wt.*cos(sample  +  P);  %  HARMONIC  CURRENT  OF  THE  LOAD 

Ic(:,k)  =  (Sh( :,  index +l).*sin(sample)  +  Ch(:,  index  +  l).*cos(sample)>; 

%  CONDITIONER  CURRENT 

Vf(k-1)  =  (If(k)  -  If(k-2))*(L*N*f/2)  +  If(k-1)*R;  %  FUNDAMENTAL  VOLTAGE 

Ve(:,k-1)  =  (Ih(:,k)  +  Ic(:,k)  -  Ih(:,k-2)  -  Ic(:,k-2)).*(L*N*f/2)  +  (Ih(:,k-1)  +  Ic(:,k-1».*R; 

%BUS  ERROR  VOLTAGE 

end 
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MSE(:,2*index  +  l)=  sqrt(mean(Ve(:,k-l-lag:k-l)’.A2))Y°MEAN  SQUARE  OF  ERROR  VOLTAGE 


thd{2*index+ 1) 
final 

GradChf:, index) 
Chf:,index+2) 


=  sqrt(sum(maxfVe(:,k-l-lag:k-l)’).*2/2)*100/(max(Vf(k-l-lag:k-l))/sqrtf2)); 

%  TOTAL  HARMONIC  DISTORTION 
=  k  +  1;  %  STEP  INDEX 

=  (MSEf:,2*index+ 1)  -  MSE(.,2*index))./(Ch(:,index+ 1)  -  Chf:, index)); 

%  COSINE  HARMONICS  GRADIENT  CALCULATION 
=  Ch(:,index  + 1)  -  GradCh(:, index). *mu. *MSE(:,2*index+ 1); 

%  PREDICTED  COSINE  WEIGHTING  FACTOR 


end 


%  %  %  %  %  %  %  %  %  %%%%  %  %  %  %  %  %  %  PLOT  SECTION  %%%%%%%%%%%%%%%%%%%%%% 

plot(Ve’);title(’Error  Voltage'); 

xlabelC Samples  (N)’);ylabel(’Voltage  (V)’);pause 

plot( Ic ’ ) ; title( ' Conditioner  Current ’ ) ; 

xlabelC  Samples  (N)’);ylabel(’Current  (I)’);pause 

plot(MSE’);grid;title(’Expectation’);grid; 

xlabelC  Samples  (N)  ’ )  ;y  label(  ’Magnitude  ’);  pause 

plot(Sh’);title(’Sin  weighting  coefficients ’);grid; 

xlabel(  ’ Samples  ( N) ’ ) ;y labelf  Magnitude ’ ) ;pause 

plot(Ch’);title(’Cosine  weighting  coefficients’ );grid; 

xlabelC  Samples  ( N)’ )  ;y  labelf  ’  Magnitude’ ) ; 

plot(GradSh’);grid;title(’Grad  Sh’); 

xlabelf  ’  Samples  (N)  ’ )  ;y labelf  ’  Magnitude  ’ )  ;pause 

plotfGradCh’);grid;titIe(’Grad  Ch’); 

xlabelf ’Samples  (N)’);ylabelf’Magnitude’); 

plot(thd);title(’Total  Harmonic  Distortion’);grid 

xlabelf’Number  of  Cycles ’);ylabelf’ Percent  (%)'); 
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APPENDIX  B  -  (MATLAB)  ADAPTIVE  NONLINEAR  MODEL 


%%  THESIS  PROGRAM  2  (NON  LINEAR  LOAD) 

%%  JOEL  ZUPFER 

%%  28  MAY  93 

%  %  REVISED  23  NOVEMBER  93 

%  %  SIMULATION  OF  CIRCUIT  FIGURE  2.4 

clear 

R  =  700;  %  MODELLED  RESISTOR  VALUE  (OHMS) 

Li  =  1.857;  %  MODELLED  INDUCTOR  VALUE  (HENERY’S)  INITIAL  (NONLINEAR  L) 

f  =  60;  %  FUNDAMENTAL  FREQUENCY  (HZ) 

N  =120;  %  NUMBER  OF  SAMPLES  PER  CYCLE 

lag  =  120;  %LAG  FOR  CALCULATION  OF  THE  MSE 

Cycle  =  input(’ Number  of  cycles  =  ’);  %  NUMBER  OF  CYCLES  IN  SIMULATION 

check  =  input(’Number  of  cycles  between  gradient  calculations  =  ’); 

%  CYCLES  BETWEEN  GRADIENT  CALCULATIONS 
totN  =  N*Cycle*check  %  TOTAL  NUMBER  OF  POINTS  IN  SIMULATION 

Delay  =  N+check;  %COUNTER  FOR  CHECKING  ERROR  VOLTAGE 

accel  =  3e-7;  %ACCELERATION  FACTOR 

Wt  =  [.51  .16  .056  .035  .025  .025  .02  .02  .015  .015]’; 

%  INITIAL  WEIGHTS  OF  ODD  HARMONICS 
H  =  [3  5  7  9  11  13  15  17  19  21]’;  %ODD  HARMONICS  OF  FUNDAMENTAL  FREQUENCY 

P  =  pi*[l/3  1/4  1/5  3/2  4/2  5/2  1/6  2/6  3/6  8/6]’;  %INITIAL  PHASE  OF  ODD  HARMONICS 

mu  =  H.~(-l).*accel;  %ACCELERATION  FACTOR  ADJUSTED  FOR  HARMONICS 


%%%%%%%%%%%%%%%% INITIALIZE  PARAMETER  MEMORY%%%%%%%%%%%%%%% 


Ih 

=  zeros(length(Wt),totN); 

%  HARMONIC  CURRENT 

Ic 

=  zeros(length(Wt),totN); 

%CONDITIONER  CURRENT 

If 

=  zeros(l.totN); 

%  FUNDAMENT  AL  CURRENT 

Ve 

=  zeros(length(Wt),totN-l); 

%BUS  VOLTAGE 

Sh 

=  zeros(length(Wt),round(Cycle/2)+2); 

%HARMONIC  SIN  WEIGHTS 

Ch 

=  zeros(length(Wt),round(Cycle/2)+2); 

%HARMONIC  COS  WEIGHTS 

MSE 

=  zeros(length(Wt),Cycle+l); 

%MEAN  SQUARE  ERROR 

thd 

=  zeros(l,Cycle+l); 

%  TOTAL  HARMONIC  DISTORTION 

GradSh 

=  zeros(length(Wt),Cycle); 

%  GRADIENT  OF  SIN  HARMONICS 

GradCh 

=  zeros(length(Wt), Cycle); 

%  GRADIENT  OF  COS  HARMONICS 

%%%%%%%%%%%%%%%%%%%  INmALCONDITIONS%  %%%%%%%%%%%%%%%%%%% 
Sh(:,2)  =  ones(length(Wt),l)./HM;  %INITIAL  WEIGHTING  FACTOR  FOR  SINE  HARMONICS 
Ch(:,2)  =  ones(length(Wt),  1)./H*.  1  ,•% INITIAL  WEIGHTING  FACTOR  FOR  COSINE  HARMONICS 
If(l)  =  10*sqrt(2;; 

If(2)  =  9.9*sqrt(2);  PEAK  CURRENT  VALUES 


34 


%%%%%%%%%%%%%%%%%% %MAIN  PROGRAM %%%%%%%%%%%%%%%%%%%%%% 

% %%%%%%%%%% %BUS  VOLTAGE  WITH  NO  CONDITIONER  CURRENT% %%%%%%%%%% 
for  k  =  3:N-1 

sample  =  H.*2*pi*k/N;  %  DISCRETE  SAMPLE  POINT 

If(k)  =  I0*sqrt(2)*cos(2*ni*k/N);  %  FUNDAMENTAL  CURRENT 

Ih(:,k)  =  Wt.*cos(sample  +  P)  %  HARMONIC  CURRENT  OF  THE  LOAD 

Ic(:,k)  =  Sh(:,l).*sin(sample)  +  Ch(:,l).*cos(sampte);  %CONDITIONER  CURRENT 

L  =  Li  *ones(10,l)  ./  (1  +abs((Ih(:,k-l)+Ic(:,k-l))/10));  %  NONLINEAR  INDUCTANCE 

Vf(k-1)  »  (If(k)  -  If(k-2))*(Li*N*f/2)  +  If(k-1)*R;  %  FUNDAMENTAL  VOLTAGE 

Ve(:,k-1)  =  (Ih(:,k)  +  Ic(:,k)  -  Ih(:,k-2)  -  Ic(:,k-2)).*(L*N*f/2)  +  (Ih(:,k-1)  +  Ic(:,k-1)).*R; 

%BUS  VOLTAGE 
end 

MSE(:,1)=  sqrt(mean(Ve(:,l:k-l)’.*2))’;  %MEAN  SQUARE  OF  ERROR  VOLTAGE 

thd(l)  =  sqrt(sum(max(Ve(:,2:k-l)’).A2/2)*100/(max(Vf(2:k-l))/sqrt(2)); 

%  TOTAL  HARMONIC  DISTORTION 
final  =  k  +  1;  %  STEP  INDEX 

%%%%%%%%%%%%  %BUS  VOLTAGE  WITH  CONDITIONER  CURRENT  %  %%%%%%%%%%% 
for  index  =  2:  Cycle 

fork  =  final:  Delay + final- 1 

sample  =  H.*2*pi*k/N;  %  DISCRETE  SAMPLE  POINT 

If(k)  =  10*sqrt(2)*cos(2*pi*k/N);  %  FUNDAMENTAL  CURRENT 

Ih(:,k)  =  Wt.*cos(sample  +  P);  %  HARMONIC  CURRENT  OF  THE  LOAD 

Ic(:,k)  =  Sh(:, index).  *sin(sample)  +  Ch(:,  index- l).*cos(sample); 

%  CONDITIONER  CURRENT 

L  =Li  *ones(10,l)  ./  (1  +  abs((Ih( : tk- 1 ) + Ic( : ,k- 1 ))/ 10)%  NONLINEAR  INDUCTANCE 

Vf(k-1)  =  (If(k)  -  If(k-2))*(Li*N*f/2)  +  If(k-1)*R;  %  FUNDAMENT  AL  VOLTAGE 

Ve(:,k-1)  =  (Ih(:,k)  +  Ic(:,k)  -  Ih(:,k-2)  -  Ic(:,k-2)).*(L*N*f/2)  +  (Ih(:,k-1)  +  Ic(:,k-1)).*R; 

%BUS  VOLTAGE 
end 

MSE( : ,  2  *index-2)  =  sqrt(mean(Ve(:,k-l-l:k-l)’/2))’; 

%MEAN  SQUARE  OF  ERROR  VOLTAGE 

thd(2*index)  =  sqrt(sum(max(  V  e( :  ,k- 1  -lag :  k- 1 )’ ) .  *2/2)*  1 00/(max(  V  f(k- 1  -lag  :k- 1  ))/sqrt(2)); 

%  TOTAL  HARMONIC  DISTORTION 
final  =  k+1;  %  STEP  INDEX 

GradSh(:, index)  =  (MSE(:,2*index-2)  -  MSE(:,2*index-3))./(Sh(:, index)  -  Sh(:, index-1)); 

%SINE  GRADIENT  CALCULATION 
Sh(:, index-1- 1)  =  Sh(:, index)  -  GradShf:, index). *mu.*MSE(:,2*index-2); 

%  PREDICTED  SINE  WEIGHTING  FACTOR 
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for  k  =  final: Delay + final- 1 

sample  =  H.*2*pi*k/N;  %  DISCRETE  SAMPLE  POINT 

If(k)  =  10*sqrt(2)*cos(2*pi*k/N);  %  FUNDAMENTAL  CURRENT 

Ih(:,k)  =  Wt.*cos(sample  +  P);  %HARM0N1C  CURRENT  OF  THE  LOAD 

Ic(:,k)  =  Sh(:, index).  *sin(  sample)  +  Ch( : ,  index).  *cos(sample); 

%  CONDITIONER  CURRENT 

L  =  Li  *ones(10,l)  ./  (l+abs((Ih(:,k-l)+Ic(:,k-l))/10)%NONLINEAR  INDUCTANCE 

Vf(k-l)  =  (If(k)  -  If(k-2))*(Li*N*f/2)  +  If(k-1)*R;  %FUNDAMENTAL  VOLTAGE 

Ve(:,k-1)  =  (Ih(:,k)  +  Ic(:,k)  -  Ih(:,k-2)  -  Ic(:,k-2)).*(L*N*f/2)  +  (Ih(:,k-1)  +  Ic(:,k-1)).*R; 

%BUS  VOLTAGE 


MSE(:,2*index-l)  =  sqrt(mean(Ve(:,k-l-l:k-l)’.^2))’; 

%MEAN  SQUARE  OF  ERROR  VOLTAGE 

thd(2*index+l)=  sqrt(sum(max(Ve(:,k-l-lag:k-l)’).A2/2)*l(X)/(max(Vf(k-l-lag:k-l))/sqrt(2)); 

%  TOTAL  HARMONIC  DISTORTION 
final  =  k+1;  %  STEP  INDEX 

GradChf index)  =  (MSE(:,2*index-l)  -  MSE(:,2*index-2))./(Ch(:, index)  -  Ch(:, index- 1)); 

%  COSINE  GRADIENT  CALCULATION 
Ch(:,index+1)=  Ch(:, index)  -  GradChf:, index). *mu.*MSE(:,2*index-l); 

%  PREDICTED  COSINE  WEIGHTING  FACTOR 


%%%%%%%%%%%%%%%% %%%PLOTTING  SECTION %%%%%%%%%%%%%%%%%%%% 

plot(Ve’);title(’Error  Voltage’); 

xlabeK’ Samples  (N)’);ylabel(’Voltage  (V)’);pause 

plot(Ic’);title(’Conditioner  Current’); 

xlabel(’Samples  (N)’);ylabel(’Current  (I)’);pause 

plot(MSE’);grid;title(’Mean  Square  Error’);grid; 

xlabelf’Samples  (N)’);ylabel(’Magnitude’);pause 

plot(Sh’);title(’Sin  weighting  coefficients’);grid; 

xlabel(  ’  Samples  (N)  ’ )  ;y label(’  Magnitude  ’ ) ; pause 

plot(Ch’);title(’Cosine  weighting  coefficients’);grid; 

xlabel(’Samples  (N)’);ylabel(’Magnitude’); 

plot(GradSh’);grid;title(’Grad  Sh’); 

xlabel(’Samples  (N)’);ylabel(’Magnitude’);pause 

plot(GradCh’);grid;title(’Grad  Ch’); 

xlabelf ’Samples  (N) ’ ) ;ylabel( ’ Magnitude ’ ) ; 

plot(thd);title(’Total  Harmonic  Distortion’);grid 

xlabeK’Number  of  Cycles’);ylabel(’Percent  (%)’); 


APPENDIX  C  -  (MATLAB)  OPTIMAL  ESTIMATION  MODEL 


%%  THESIS  PROGRAM  3 
%%  JOEL  ZUPFER 
%%  1  JUNE  93 

%  %  REVISED  30  NOVEMBER  93 
%  %  SIMULATION  OF  CIRCUIT  FIGURE  2.3 


clear 

R 

L 

f 

Wt 

H 

P 

M 

ALPHA 

Cycle 


=  700;  %  MODELED  RESISTOR  VALUE  (OHMS) 

=  1.857;  %  MODELED  INDUCTOR  VALUE  (HENERY’S) 

=  60;  %  FUNDAMENTAL  FREQUENCY  (HZ) 

=  10  0  .51  0  .16  0  .056  0  .035  0  .025  0  .025  0  .02  0  .02  0  .015  0  .015]; 

%  WEIGHTS  OF  THE  HARMONICS  INITIAL  CONDITION 
=  ]1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21]; 

%  HARMONICS  OF  FUNDAMENTAL  FREQUENCY 
=  pi*[0  0  1/3  0  1/4  0  1/5  0  3/2  0  4/2  0  5/2  0  1/6  0  2/6  0  3/6  0  8/6]; 

%  PHASE  OF  HARMONICS  INITIAL  CONDITION 
=  2*length(H);  %  NUMBER  OF  SAMPLE  POINTS/PERIOD 

=  0.75;  %  WEIGHTING  FACTOR  FOR  CONTROL  DETERMINATION 

=  input(’ Number  of  cycles  =  ’);  %  NUMBER  OF  CYCLES  IN  SIMULATION 


%%%%%%%%%%%%%%  %  INITIALIZE  PARAMETER  MEMORY  %  %%%%%%%%%%%%% 


D 

HL 

Qhat 

Q 

thetaO 

theta42 

theta 

GO 

G42 

G 

If 

Vf 

Ve 

X 

E 

Ek 

Ehat 

Uk 

Uhat 

Yk 

thd 


=  [1:42];  %  DIAGONAL  ELEMENTS  OF  H 

=  diag(D);  %  TRANSFER  MATRIX  OF  SYSTEM  LOAD 

=  eye(M)M ;  %INITIAL  GUESS  FOR  INVERSE  OF  H 

=  zeros(2*length(H),2*length(H));  %INITIAL  VALUE  OF  Q  MATRIX 

=  zeros(l.l);  %THE  ESTIMATION  FOR  THE  DC  COMPONENT 

=  zerosd.l);  %THE  ESTIMATION  FOR  THE  21st  HARMONIC 

=  zeros(2,2*(length(H)- 1 ));  %THE  ESTIMATION  OF  THE  SYSTEM  IMPEDANCE 

=  50;  %  INITIAL  VALUE  OF  GAIN  MATRIX(MEASURE  OF  UNCERTAINTY) 

=  50;  %  INITIAL  VALUE  OF  GAIN  MATRIX 


=  [ones(2*H(length(H)-l),l),zeros(2*H(length(H)-l),2),ones(2*H(length(H)-l),l)]*50; 


=  zeros(l.M); 

=  zeros(l.M); 

=  zeros(l.M); 

=  zeros(M.M); 

=  zeros(M, Cycle); 
=  zeros(M, Cycle); 
=  zeros(M, Cycle); 
=  zeros(M.Cycle); 
=  zeros(M,  Cycle); 
=  zeros(M.Cycle); 
=  zeros!  1, Cycle); 


%  INITIAL  VALUE  OF  GAIN  MATRIX 
%  INITIALIZE  THE  FUNDAMENTAL  CURRENT 
%  INITIALIZE  THE  FUNDAMENTAL  VOLTAGE 
%  INITIALIZE  THE  NOISE  VECTOR 
%  INITIALIZE  THE  MATRIX  OF  SUM  OF  COS  &  SIN 

%  ACTUAL  ERROR 
%  PERIODIC  ERROR 
%  ESTIMATED  PERIODIC  ERROR 
%  PERIODIC  CONTROL  INPUT 
%  ESTIMATED  PERIODIC  CONTROL  INPUT 
%  PERIODIC  SYSTEM  OUTPUT 
%TOTAL  HARMONIC  DISTORTION 


37 


%%%%%%%%%%%%%%%%%%% %MAIN  PROGRAM% % %%%%%%%%%%%%%%%%% 
for  n  =  0:M-1 

If(n+1)  =  10*sqrt(2)*cos(2*pi*n); 

Sample  =  H.*(2*pi*n/M); 

X(:,n  +  1)  =  reshape<[cos(Sample);sin(Sample)J,2*length(H),l); 

%  MATRIX  FOR  SAMPLED  VALUE  AS  SUM  OF  COS  &  SIN 
end 


W  =  reshape([cos(P).*Wt;sin(P).*Wt],2*length(H),l); 

%  HARMONIC  CURRENT  WEIGHT  VECTOR 
Ih  =  W’*X;  %THE  HARMONIC  CURRENT  OF  ONE  PERIOD 

Vf(l)  =  (If(2)  -  If(M))*(L*M*f/2)  +  If(l)*R;  %  FUNDAMENTAL  VOLTAGE 

Ve(:,l)  =  (Ih(:,2)  -  Ih(:,M)).*(L*M*f/2)  +  ih(:,l).*R;  %  ERROR  VOLTAGE 


for  n  =  3:M 

Vf(n-I)  =  (If(n)  -lf(n-2))*(L*M*f/2)  +  !f(n-l)*R;  %  FUNDAMENTAL  VOLTAGE 

Ve(:,n-I)  =  (Di(:,n)  -  ni(:.n-2)).*(L*M*f/2)  +  Ih(:,n-1).*R;  %  ERROR  VOLTAGE 

end 


Vf(M) 

Ve(:,M) 

MAGf 

Yk(:.l) 

E(:,l) 

Ek(:.l) 

Uk(:,2) 


=  (If(l)  -  If(M-l))*(L*M*f/2)  +  lf(M)*R;  %  FUNDAMENTAL  VOLTAGE 

=  (Ih(:,l)  -  Ih(:,M-l)).*(L*M*f/2)  +  Ih(:,M).*R;  %  ERROR  VOLTAGE 

=  max(V0/SQRT(2);  %RMS  VALUE  OF  THE  FUNDAMENTAL 

=  H*Uk(:,l);  %  SYSTEM  OUTPUT 

=  Ve’;  %  ACTUAL  ERROR 

=  (fftrig(E(:,l)))’;  %PERIODIC  ERROR  COEFFICIENTS  FREQUENCY  DOMAIN 

=  Uk(:,l)  -  ALPHA*Qhat*Ek(:,l);  %  CONTROL  COEFFICIENT  UPDATE 


fori  =  l:length(Ek(:,l))/2  -1 

MAG(i,l)  =  sqrt(sum(Ek(2*i:2*i + 1 , 1)). A2»;  %SUM  OF  THE  ERROR  VOLTAGES 

end 


MAG(length(Ek(:,l))/2,l)  =  Ek(length(Ek(:,l)),D; 

thd(l)  =  sqrt(sum(MAG.*2/2))*100/(MAGD;  %  TOTAL  HARMONIC  DISTORTION 


for  k  =  2:Cycle 

Ek(:,k)  =  fftrig(Ve)’  +  H*Uk(:,k); 

E(:,k)  =  ifftrig(Ek(:,k»’; 

Ehat(:,k)  =  Ek(:,k)  -  Ek(:,k-1); 
Uhat(:,k)=  Uk(:,k)  -  Uk(:,k-1); 


%  RECURSIVE  LOOP  FOR  ERROR  CALCULATION 
%ACTUAL  ERROR  IN  FREQUENCY  DOMAIN 
%TIME  DOMAIN  OF  ERROR 
%  ESTIMATED  ERROR 
%  ESTIMATED  CONTROL 


thetaO 

GO 


=  thetaO  +  (Uhat(l.k)  -  Ehat(  1  ,k)*thetaO)*GO*Ehat(  1  ,k)/(  1  +  Ehat(l,k)*GO*Ehat(l,k)); 

%  RECURSIVE  LEAST  SQUARES  ESTIMATION  OF  COEFFICIENT 
=  GO  -  GO*Ehat(l,k)*Ehat(l,k)*GO/(l  +  Ehat(l,k)*GO*Ehat(l,k)); 

%  UPDATE  OF  THE  GAIN 


Q(  1,1)  =  thetaO; 
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for  h  =  l:2*H(length(H)-l) 

gain  =  reshape(G(h,:),2,2);  %GAIN  MATRIX  IS  2X2 

if  rem(h  ,2)  =  =  l;  %  FOR  THE  ODD  LOOPS 

Phi  =  Ehat(h+  l:h+2,k)’.*[l  -1|; 
else 

Phi  =  flipud(Ehat(h:lHT,k))’;  %FOR  THE  EVEN  LOOPS 

end 

theta(:,h)  =  theta(:,h)  +  (Uhat(h+l,k)  -  Phi*theta(:,h))*gain*Phi7(l  +  Phi*gain*Phi'); 

%  RECURSIVE  LEAST  SQUARES  ESTIMATION  OF  COEFFICIENTS 
gain  =  gain  -  gain*Phi’*Phi*gain/(l  +  Phi*gain*Phi’);  %  UPDATE  OF  THE  GAIN 

G(h,:)  =  reshape(gain,l,4);  %CHANGE  GAIN  MATRIX  BACK  TO  1X4 

Q(h+l,h+l:h+2)  =  theta(:,h)';  %  PLACE  ESTIMATES  WITH  IN  THE  Q  MATRIX 


end 


theta42  =theta42  +  (Uhat(42,k)  -  Ehat(42,k)*theta42)*G0*Ehat(42,k)/... 

(I +Ehat(42,k)*G0*Ehat(42,k)); 

%  RECURSIVE  LEAST  SQUARES  ESTIMATION  OF  COEFFICIENT 
G42  =  G42  -  G42*Ehat(42,k)*Ehat(42,k)*G42/(l  +  Ehat(42,k)*G42*Ehat(42,k)); 

%  UPDATE  OF  THE  GAIN 


Q(42,42) 
Uk(:,k+  L) 
Yk(:,k) 
Nfk 


theta42;  %PLACE  ESTIAMTES  WITH  IN  THE  Q  MATRIX 

Uk(:,k)  -  ALPHA*Q*Ek(:,k);  %  UPDATE  THE  CONTROL 

ifftrig(H*Uk(:,k))’;  %  SYSTEM  OUTPUT 

fftrig(Ve  +  Yk(:,k)’)’;  %  ERROR  AFTER  NEW  CONTROL  INPUT 


for  i  =  1  :length(Nfk(:,k))/2  -  1 

MAG(i.k)  =  sqrt(sum(Nfk(2*i:2*i  +  l),k)."2));  %  MAGNITUDE  OF  ERROR  VECTOR 

end 


MAG(length(Nfk(: ,k))/2,k)  =  Nfk(length(Nfk(:,k)),k);  %SUM  OF  THE  ERROR  VOLTAGES 
thd(k)  =  sqrt(sum(MAG(:,k).A2/2))*100/<MAGf);  %  TOTAL  HARMONIC  DISTORTION 


%%%%%%%%%%%%%%%%% % %PLOT  SECTION %%%%%%%%%%%%%%%%%%%%% 

plot(E’);title(’Actual  Error’);gnd; 

xlabel(’Samples  (N)’);ylabel(’Voltage  (V)’);pause 

plot(Ek’);title(’ Periodic  Error’);grid; 

xlabel(’ Samples  (N)’);ylabel(’Voltage  (V)’);pause 

plot(Ehat’);title(’Estimated  Error’);grid 

xlabel(’Samples  (N)’);ylabel(’Voltage  (V)’);pause 

plot(Uk’);title(’Control  weighting  coefficients’);grid; 

xlabelf’Samples  (N)’);ylabel(’Magnitude’);pause 

plot(Uhat’);title(’Estimated  Control  weighting  coefficients’);grid; 

xlabelf  Samples  (N)’);ylabel(’Magnitude’);pause 

plot(Yk’);title(’System  output’);grid; 

xlabelf ’Samples  (N)’);ylabel(’Magnitude’>; 

plot(thd);  titlef Total  Harmonic  Distortion’);grid 

xlabelf’Number  of  Periods  (k)’);ylabel(’Percent  (%)’); 
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(function  w  =  fftrig(x) 

%  w=fftrig(x) 

%  it  computes  the  fft  coefficients  of  the  vector  x 
%  in  trigonometric  form. 

%  x(n) = w  l  +  w2*cos(2pi/N) 4- w3*sin(2pi/N) + w4*cos(4pi/N) + w5*sin(4pi/N)  + . . . 
%  ...  +  w(N-2)cos(2pi(N/2- 1  )/N) + w(N- l)sin(2pi(N/2- 1 )/N) + w(N)cos(pi  n); 

%  where  N =length(x),  assumed  to  be  a  power  of  2  (if  not  x  is  padded  with  0’s) 

X  =  fft(x); 

N  =length(X); 

Xr  =real(X); 

Xi  =imag(X); 
w(  1)  =Xr(l)/N; 
k=l:l:(N/2)-l; 

w(2*k)  =2*Xr(2:(N/2»/N; 

w(2*k+l)  =-2*Xi(2:(N/2))/N; 
w(N)  =Xr((N/2)  +  l)/N; 

end  %  fftrig 


function  x=ifftrig(w) 

%  x=iffirig(w) 

%  compute  time  sequence  x  from  trig,  coefficients  w. 
%  See  FFTRIG 
N  =length(w); 

Xr(l)  =w(l)*N; 

Xi(l)  =0; 
k=l:l:(N/2)-l; 

Xr(2:(N/2))  =N*w(2*k)/2; 

Xi(2:(N/2))  =-N*w(2*k+l)/2; 

Xr((N/2)+l)  =N*w(N); 

Xi((N/2)+l)  =0; 

X  =Xr+sqrt(-l)*Xi; 

X(N-k+l)  =conj(X(k+l)); 
x  =real(ifft(X)); 

end  %  ifftrig 
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